/*
 * pulseseq.cpp
 *
 *  Created on: May 26, 2011
 *      Author: nzhao
 */

#include "share/pulseseq.h"

int		CDD_ORDER;
int     SEQ_NAME;

double	* PULSE_PACK_POSITION;
int		MAX_PULSE_IN_PACK = 100;
int		PACK_PULSE_NUM;

int		INTERVAL_NUM;
double	* TIME_SEQ_INTERVAL;
double	* PULSE_TIMING;
int		* PULSE_SYMMETRY_MODE;

int set_time_sequence()
{
    // This function set the value of global variable TIME_SEQ_INTERVAL,
    // which is a INTERVAL_NUM * (NTIME+1) array.

    int		i, j, count;
    double	timeval;

    //////////////////////////////////////////////////////////////////
    //determine the pulse sequence
    if (strcmp(STR_SEQ_NAME, "UDD") == 0) //parameter 1 == UDD
    {
    	SEQ_NAME = 1;
    }
    else if (strcmp(STR_SEQ_NAME, "PDD") == 0) //parameter 2 == PDD
    {
    	SEQ_NAME = 2;
    }
    else if (strcmp(STR_SEQ_NAME, "CDD") == 0) //parameter 3 == CDD
    {
    	SEQ_NAME = 3;
    }
    else if (strcmp(STR_SEQ_NAME, "Pack") == 0) //parameter 4 == Pack
    {
    	SEQ_NAME = 4;
    }
    else
    {
    	SEQ_NAME = -1;
    	cout << "***** ERROR: unknown pulse sequence. *****" << endl;
    	assert(0);
    }

    //////////////////////////////////////////////////////////////////
    // Special treatment for CDD and packaged pulse
    if (SEQ_NAME == 3)   //CDD
    {
    	CDD_ORDER = PULSE_NUM;
    	PULSE_NUM = get_CDD_pulse_num(CDD_ORDER);
    }
    if (SEQ_NAME == 4)	//Pack
    {
    	import_pulse_pack();
    	PULSE_NUM = PULSE_NUM * PACK_PULSE_NUM;
    }

    //////////////////////////////////////////////////////////////////
    // Calculate pulse intervals
    // This function set the value of global variable TIME_SEQ_INTERVAL,
    // which is a INTERVAL_NUM * (NTIME+1) array.

    INTERVAL_NUM = 2 * ( PULSE_NUM + 1);

    TIME_SEQ_INTERVAL	= new double[ INTERVAL_NUM * (NTIME+1)];
    PULSE_TIMING		= new double[ PULSE_NUM + 2];

    count = 0;
    for (i = 0; i <= NTIME; i++)
    {
            timeval = MIN_TIME + i * TIME_STEP;

            PULSE_TIMING[0] = 0.;
            for (j = 1; j <= PULSE_NUM; j++)
            {
                    PULSE_TIMING[j] =  getPulse(j);
            }
            PULSE_TIMING[PULSE_NUM + 1] = 1.0;


            for (j = 0; j < PULSE_NUM + 1; j++)
            {
                    TIME_SEQ_INTERVAL [i * INTERVAL_NUM + j] = timeval * ( PULSE_TIMING[j + 1] - PULSE_TIMING[j] );
                    count ++;
            }
            for (j = 0; j < PULSE_NUM + 1; j++)
            {
                    TIME_SEQ_INTERVAL [i * INTERVAL_NUM + PULSE_NUM + 1 + j] = TIME_SEQ_INTERVAL [i * INTERVAL_NUM + PULSE_NUM - j ];
                    count ++;
            }
    }
    assert( count == INTERVAL_NUM * (NTIME+1) );

//////////////////////////////////////
// set pulse symmetry mode
// PULSE_SYMMETRY_MODE[0] >   0: PDD
// PULSE_SYMMETRY_MODE[0] ==  0: symmetric pulse, i.e. t_i = t - t_(N-i+1)
// PULSE_SYMMETRY_MODE[0] == -1: asymmetric pulse.
    int quot, md, bitlen = 0, * temp, np;
    temp = new int [20];
    PULSE_SYMMETRY_MODE = new int [20];

    switch(SEQ_NAME)
    {
    case 1://UDD
  		PULSE_SYMMETRY_MODE[0] = -1;
    	break;
    case 2://PDD
    	if(PULSE_NUM < 5)
    	{
    		PULSE_SYMMETRY_MODE[0] = -1;
    	}
    	else
    	{
    		quot = PULSE_NUM / 2;
    		do
    		{
    			md = quot % 2;
    			temp[bitlen] = md;
    			bitlen ++;
    			quot /= 2;
    		}
    		while(quot != 0);

    		for(i = 1; i <= bitlen; i ++)
    			PULSE_SYMMETRY_MODE[i] = temp[i-1];
    		PULSE_SYMMETRY_MODE[0] = bitlen;

    		cout << "PDD Pulse Symmetry:" << PULSE_SYMMETRY_MODE[0] << endl;
    		for(i = 1; i <= bitlen; i++)
    			cout << PULSE_SYMMETRY_MODE[i] << endl;
    		cout << endl;
    	}
    	break;
    case 3://CDD
    	PULSE_SYMMETRY_MODE[0] = -1;
    	break;
    case 4: //Packaged DD
    	PULSE_SYMMETRY_MODE[0] = -1;
    	break;
    default:
    	PULSE_SYMMETRY_MODE[0] = -1;
    }

    return 1;
}

double getPulse(int n)
{
	double res;

	switch(SEQ_NAME)
	{
	case 1:
		res = sin( n * PI / 2. / ( PULSE_NUM + 1 ) );
		res *= res;
		break;
	case 2:
		res = ( 2.0 * n - 1.) /( 2. * PULSE_NUM );
		break;
	case 3:
		res = cddpulse(CDD_ORDER, n);
        break;
	case 4:
		res = packpulse(n);
		break;
	default:
		res =  1.;
	}
	return res;
}

double packpulse(int n)
{
	double res, pack_num;

	int pack_i, idx_inpack;

	pack_i = (n - 1) / PACK_PULSE_NUM;
	idx_inpack = ( n - 1 ) % PACK_PULSE_NUM;

	pack_num = (double) (PULSE_NUM / PACK_PULSE_NUM);
	res = ( (double)pack_i + PULSE_PACK_POSITION[idx_inpack] ) / pack_num;

	return res;
}

double cddpulse(int cdd_order, int k)
{
    double res, lower_order_pos;
    int npulse, half_n, q;

    npulse = get_CDD_pulse_num(cdd_order);

    if( cdd_order == 1)
    {
        res = 0.5 * k;
        return res;
    }
    if( npulse % 2 == 1 && k == (npulse + 1 ) / 2 )
        return 0.5;

    half_n = npulse % 2 == 0 ? npulse / 2 : ( npulse + 1) / 2;
    q = k > half_n ? k - half_n : k;

    lower_order_pos = 0.5 * cddpulse(cdd_order - 1, q);

    if( cdd_order % 2 == 0)
    {

        if( k <= npulse /  2 )
        {
            res = lower_order_pos;
        }
        else
        {
            res = 0.5 + lower_order_pos;
        }
    }
    else
    {
        if( k <= (npulse - 1) /2 )
        {
            res = lower_order_pos;
        }
        else if ( k == (npulse + 1) / 2 )
        {
            res = 0.5;
        }
        else
        {
            res = 0.5 + lower_order_pos;
        }
    }
    return res;
}

int get_CDD_pulse_num(int cdd_order)
{
    if(cdd_order == 1)
        return 1;
    if(cdd_order%2 == 0)
    {
        return 2 * get_CDD_pulse_num(cdd_order-1);
    }
    else
    {
        return 2 * get_CDD_pulse_num(cdd_order-1) + 1;
    }

}

void import_pulse_pack()
{
	int npulse, i;
	char packfilename [300];

	strcpy(packfilename, PULSE_PATH);
	strcat(packfilename, INPUT_PACK_NAME);
	strcat(packfilename, ".ppk");
	cout << "Pack file: " << packfilename << endl;
	ifstream pulse_pack(packfilename);
	if(!pulse_pack)
	{
		cout << "Cannot open pulse file." << endl;
		assert(0);
	}

	pulse_pack >> npulse;
	assert(npulse < MAX_PULSE_IN_PACK);

	PACK_PULSE_NUM = npulse;
	PULSE_PACK_POSITION = new double [npulse];

	for(i = 0; i < npulse; i++)
	{
		pulse_pack >> PULSE_PACK_POSITION[i];
	}

	pulse_pack.close();
}
